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ABSTRACT 

Torque fluctuations due to magnetorotational turbulence in proto-planetary disks 
may greatly influence the migration patterns and survival probabilities of nascent plan- 
ets. Provided that the turbulence is a stationary stochastic process with finite amplitude 
and correlation time, the resulting diff'usive migration can be described with a Fokker- 
Planck equation, which we reduce to an advection-diff'usion equation. We calibrate 
the coefficients with existing turbulent-disk simulations and mean- migration estimates, 
and solve the equation both analytically and numerically. Diffusion tends to dominate 
over advection for planets of low-mass and those in the outer regions of proto-planetary 
disks, whether they are described by the Minimum Mass Solar Nebula (MMSN) or by 
T-Tauri alpha disks. Diffusion systematically reduces the lifetime of most planets, yet it 
allows a declining fraction of them to survive for extended periods of time at large radii. 
Mean planet lifetimes can even be formally infinite (e.g. in an infinite steady MMSN), 
though median lifetimes are always finite. Surviving planets may linger near specific 
radii where the combined effects of advection and diffusion are minimized, or at large 
radii, depending on model specifics. The stochastic nature of migration in turbulent 
disks challenges deterministic planet formation scenarios and suggests instead that a 
wide variety of planetary outcomes are possible from similar initial conditions. This 
would contribute to the diversity of (extrasolar) planetary systems. 

Subject headings: accretion, accretion disks — planetary systems: formation — plane- 
tary systems: proto-planetary disks 



1. Introduction 

A decade of extrasolar planet discoveries has shown that the process of planet formation is more 
complex than originally anticipated. It leads to a remarkable diversity of planetary configurations, 
ranging from migrating hot Jupiters to eccentric giant planets, as well as our own "circular" Solar 



- 2 - 



System. Following these discoveries, progress has been made in understanding planet formation, 
but the theory is still incomplete (e.g. Marcy et al. 2000; Wuchterl et al. 2000; Papaloizou &; 
Terquem 2006; Armitage & Rice 2006). 

The leading scenario for the formation of giant planets is the core accretion mechanism. Icy 
planetesimals located beyond the snow line in their host disk collide repeatedly to grow a core with 
a modest gaseous atmosphere. If this core succeeds in reaching a critical mass of a few tens of 
Earth masses, runaway accretion of a massive gaseous envelope proceeds and leads, ultimately, to 
the formation of a gaseous giant planet (Safronov 1969; Mizuno 1980; Pollack et al. 1996). Some 
evidence supporting this scenario has emerged in recent years, in the form of a metallicity trend 
for stars hosting planets (e.g. Santos et al. 2004; Gilliland et al. 2000; Weldrake et al. 2005), the 
discovery of a high density hot Jupiter (Sato et al. 2005) and that of a surprisingly low-mass micro- 
lensing planet (Beaulieu et al. 2006). However, a long-standing difficulty for the core accretion 
scenario, which has not yet been fully elucidated, is the fact that the timescale required to build-up 
critical core masses and thus large gaseous envelopes (~ 10^-10^ yr) is comparable to the lifetimes 
of proto-planetary disks (e.g. Pollack et al. 1996; Papaloizou & Terquem 1999). 

Orbital migration adds a layer of complication to theories of planet formation. As a result of 
gravitational interactions with their gaseous disk (Goldreich & Tremaine 1980; Lin &; Papaloizou 

1986), the orbits of planets in the terrestrial mass range are predicted to decay on timescales 
[r^ 10^ yr) short compared to disk lifetimes (e.g. Korycansky & Pollack 1993; Ward 1997a, b). 
Migration is slower for planets of much smaller or much larger masses: in the first case because the 
torque causing migration is quadratic in planet mass, and in the second case because the planet 
opens a gap and then migrates on the disk's accretion timescale, which can be comparable to its 
lifetime. While it is possible or probable that many terrestrial planets form by agglomeration of 
smaller bodies after the gas is gone, this is not an option for the solid cores of Jovian planets since, 
in the core accretion scenario, the cores must form before the gaseous envelopes can be accreted. 
The prevalence of Jovian planets with orbital periods of only a few days deepens the mystery as 
it suggests that these planets did migrate but stopped short of merging with their stars at orbital 
radii where even the accretion timescale would seem to have been very short (e.g. Lin et al. 1996). 

Analytic calculations and most hydro-dynamical simulations of migration usually assume a 
disk that is laminar apart from the waves and shocks excited by the planet itself (e.g. Korycansky 
& Pollack 1993; Ward 1997a; Kley et al. 2001; Nelson & Benz 2003a,b; D'Angelo et al. 2003; Lufkin 

et al. 2004; Schaefer et al. 2004). But the effective viscosity of disks probably involves turbulence. 
Nelson & Papaloizou (2004), Laughlin et al. (2004), and Nelson (2005, hereafter N2005) have 
found in 3D simulations of magneto-rotational turbulence that the instantaneous torque exerted 
on a planet in the terrestrial mass range is subject to fluctuations many times its mean value, 
apparently caused by turbulent density fluctuations in the planet's vicinity. In fact, no obvious 
secular decay manifests itself in the orbits of planets with mass Mp < lOM®, although because the 
simulations are limited to ~ 10^ planetary orbits — compare > 10^ for Jupiter during the lifetime 
of the proto-solar nebula — and the predicted decay in semimajor axis is < 10% over this period. 
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the expected secular trend might be difficult to discern amid the fluctuations. 
The possibility that planetary torques are stochastic raises several issues: 

1. Is the time-averaged torque exerted on a planet of a given angular momentum affected? N2005 
finds that turbulence excites planetary eccentricity, e. In linear theory, eccentricity allows the 
planet to interact with the disk at additional resonances that do not couple to circular orbits 
(Goldreich & Tremaine 1980; Goldreich & Sari 2003). Also, at second and higher orders in 
e, the radial positions of the original resonances shift if one compares eccentric and circular 
orbits at fixed angular momentum; they do not shift at fixed orbital energy since keplerian 
energies and mean motions vary only with semi-major axis. Thus for eccentric planets, one 
must distinguish the rate of change of the planetary angular momentum (torque) from the 
rate of change of orbital energy. Smoothing of the planetary potential along epicycles tends 
to reduce the strength of the highest-m Lindblad resonances that dominate the torque on 
circular orbits. All of these perturbative effects are expected to become significant when 
e > h/r, where h = Cg/^ is the disk thickness. Planets of terrestrial mass achieved such 
eccentricities in N2005's simulations. On the other hand, at about the same eccentricity the 
epicyclic velocities become supersonic, which will probably raise significant shocks, a non- 
perturbative effect. Thus both the sign and magnitude of the changes in the mean torque 
due to turbulently excited eccentricity are unclear. We are not aware of any systematic 
investigation of this issue, nor have we attempted one. Throughout this paper, we treat the 
planetary orbit as circular. 

2. Even for circular orbits, stochastic changes in semi-major axis on timescales long compared 
to the period but short compared to the nominal migration time may significantly alter the 
mean orbital lifetime. This is particularly likely when the local mean torque varies strongly 
with radius, as it does in the alpha disks studied by Menou & Goodman (2004). Because of 
rapid changes in opacity in certain temperature regimes, the disks exhibit sharp peaks in the 
local migration timescale. Under laminar conditions, orbital lifetimes can be dominated by 
the time required to drift across these peaks. Turbulent fluctuations in the torque, however, 
could allow planets to "jump over" the peaks, thus shortening the lifetime. On the other 
hand, diffusion to large radii where the local migration time is long might prolong the orbital 
lifetime. 

3. Finally, orbital diffusion should allow the occasional planet to survive much longer than 
the mean lifetime defined by its birthplace. If sufficiently many planets are formed, this 
might in principle help to reconcile theoretical predictions of rapid mean migration with the 
observations at least of the longer-period planets (though it seems unlikely to explain hot 
Jupiters) . 

The methods adopted in this paper allow us to explore the second and third issues above, 
but not the first. In §2, we develop a simple advection-diffusion or Fokker-Planck equation for the 
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probability distribution function of planetary orbital angular momentum in the presence of torque 
fluctuations. Paramctrizations of the advection and diffusion coefficients are calibrated against 
semianalytic formalisms and the simulations of N2005, respectively. §3 explores analytic steady- 
state and time-dependent solutions to this equation in steady disks whose properties are power 
laws in radius. §4 and §5 present numerical, time-dependent solutions in such a disk (the Minimum 
Mass Solar Nebula) and in alpha disks with realistic opacities. A summary of our main results and 
conclusions is given in §6. 



2. Advection-diffusion Equation for Planetary Migration 

Planets migrate as a result of torques exerted by the disk. Because of turbulence, the torque 
is a stochastic function with a mean part T and a fluctuating part 5T. For a given planetary mass 
(Mp) and a given time-averaged surface density and thickness of the disk [S(r), h{r) respectively], 
the mean torque depends only on the semimajor axis, r^, of the planetary orbit, which is assumed 
to be circular. The diffusion equation is most conveniently derived in terms of the orbital angular 
momentum^ J = Mpy^GM*/^, rather than Vp, as the spatial independent variable. 

The fluctuating part of the torque depends on time as well as J: ST = ST{t, J). By construction, 
Sr = 0. The fluctuating torque is taken to be a temporally stationary stochastic process with a 
finite amplitude and correlation time, so that the integral 

oo 

D{J) = ^ J dr{t - r/2, J)ST{t + r/2, J) dr (1) 

— oo 

exists. D{J) will play the role of diffusion coefficient. The correlation time can then be defined by 

T,^ D{J)/[5r{t,J)f. (2) 

Temporal stationarity implies that D{J), [Sr{t, J)]^, and Tc are independent of t. 

The desired advection-diffusion equation derives from a Fokker-Planck equation, and its formal 
derivation is completely standard (e.g. van Kampen 1981), but it is useful to repeat the derivation 
here in order to emphasize the underlying assumptions and to ask whether they are justified in this 
case. Two important assumptions have already been introduced: the finiteness of Eqs. (1) & (2). 

Let P{AJ \ At, J) be the probability that a planet of initial angular momentum J suffers a 
change A J during time interval At. Moments of the change are 

oo 

(AJ)^= J (AJ)"P(AJ|Ai, J)dAJ. (3) 

— oo 



^Throughout our analysis, we focus on the hmit Mp <C M*, where M, is the mass of the central star. 
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Of course, except for (AJ)'^ = 1, these depend implicitly on J and Ai. 



Next, let f{t, J) be the probability that the planet has angular momentum J at t > given 
some initial condition at t = 0. The probability at t + Af is 



oc 

f{t + At,J) = J P{AJ\At,J - AJ)f{t,J - AJ)dAJ 



-OD 
OO 



^ i-^(AJ)«^ [P(AJ |Ai, J) fit, J)] dAJ 



— OO 
OO 

E 

n=0 



n=0 

(-1)" a" 



n 



{(AJ)^/(i,J)} 



(4) 



In the Fokker-Planck approximation, the changes are presumed small compared to J itself but 
cumulative, so that one can ignore the terms n > 2 in the sum above. More will be said about the 
justification for this below. 



Clearly {AJf = 1, AJ = V{J)At, and 



[AJf 



At 



T{J)At + J Sr{t')dt' 



At 



'2tr} 



= TiJfiAtf + j dtm j dT 5T{tm- T/2,J)5T{tm + t/2,J). 

-2t^ 

The next important assumption is that is short compared to the timescale on which J changes 
by order itself — let us write tj for the latter timescale — so that wc may take At in the range 
Tc <C At <C tj. In that case, the double integral above ~ 2D{J)At. At this point. 



fit + At, J) ^ fit, J) + 



At + O [iAtf] , 



where the 0[(Ai)^] part includes the (TAt)^ contribution to (AJ)^. Neglecting these terms on the 
grounds that At <^tj, one has the desired advection-diffusion equation: 



df dFj 

dt dj 







Fj = f(J)/(t,J) - 



d_ 
dJ 



DiJ)fit,J) 



(5) 



An interesting feature of the above result is that if the diffusion coefficient -D( J) varies with 
J, then it influences the mean migration rate. That is to say, one can rewrite the flux in Eq. (5) as 



/- dD\ df 
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so that —dD/dJ contributes to the advection term. 



2.1. 



Calibration of the diffusion and advection coefficients 



We have tried to abstract a value and scahng for D{J) from the MHD simulations of N2005; 
the results of Laughlin et al. (2004) were difficult to translate into our framework. Nelson reports 
that the typical timescale of fluctuations in 5T is of order half an orbital period. Hence we take 



A direct estimate of the correlation function of the data suggests, however, that correlations 
may persist to much longer than an orbital period, and it is conceivable that the integral (1) for 
D{J) does not even exist (Nelson 2005, private communication). If the correlation time does not 
exist, then the methods of this paper are inapplicable to migration. Such long-term correlations 
imply a "memory" in the turbulence, presumably involving persistent structures in addition to 
the planet of interest: for example, long-lived vortices (e.g. Promang & Nelson 2005; Johansen & 
Klahr 2005), or of course other planets. Pending further numerical evidence, we adopt as a working 
hypothesis that does exist but with uncertain magnitude. 

To complete the diffusion coefficient, we require a paramctrization of the variance of the fluc- 
tuating torque in terms of time-averaged disk properties. This entails some guesswork, as the 
simulations of N2005 and Nelson &: Papaloizou (2004) have explored only a limited range of disk 
models and planetary radii. A natural scale for the gravitational force exerted on the planet by 
the local gas is 27rGEMp, where S is the surface density of the disk. This is the force that the 
planet would feel if suspended just above the disk. The corresponding natural scale for the torque 
is 27rGSMprp. Therefore, we postulate that the variance of the torque fluctuations is 



where Cd is a dimensionless coefficient depending upon the strength of the turbulence. 

Since the parametrization (7) underlies all of our results, it is worth more discussion. The tur- 
bulence is accompanied by density perturbations of r.m.s. amplitude 5p\ on scale A. An individual 
perturbation of this scale exerts a gravitational force per unit mass ~ G5p\\^d^^ at distance d 
ii X < h, and ~ G6pxhX^d~^ ii X > h. Presuming that the forces from different perturbations 
add in quadrature, the mean-square force is dominated by the closest such perturbation, at d ~ A, 
since the number of A-scale perturbations within distance d increases more slowly with d than 
the mean-square contribution from individual perturbations decreases: oc {d/X)^ for d < h and 
oc {d/Xy for d > h, versus d~^. Thus, the r.m.s. force due to scale A is Sf\ ~ G5p\X for A < /i and 
<^/a ^ Gdpxh ^ G5T,x for X> h. The contributions to D{J) from X < h are probably unimportant 




(6) 



[ST{t,J)]^ ^ {Cd X 2TrGErpMpf 



(7) 
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because of the weighting by A and because smaller scales are likely to have shorter correlation times. 
If the density power spectrum were as "blue" as white noise, 5px oc A~^/^, and if the correlation 
time Tc^A Of A, then there would be equal contributions from all scales < /i, with the result that the 
total contribution from these scales would be larger than the contribution from A > /i by a loga- 
rithmic factor. In reality, the density power spectrum is surely red, since local magnetorotational 
(henceforth MRI) simulations find that most of the turbulent magnetic and kinetic energy resides 
at scales > h (Hawley, Gammie, & Balbus 1995; Brandenburg ct al. 1995). Thus, the dominant 
contribution to the force comes from scales > h. Since 6f\ ~ GST^x in this regime, the total r.m.s. 
force fluctuation is proportional to the total r.m.s. surface-density fluctuation, 5E. 

It remains to discuss how 5p/ p and hence ^S/S may depend upon the strength of the turbu- 
lence. Density fluctuations may arise from fluctuations in gas pressure or in entropy. Local MRI 
simulations indicate that 6p/ p ~ 6P/P ~ dB^/STrpc^ « 2a (Sano et al. 2003), as might be expected 
from equipartition, since the internal energy per unit mass associated with a density fluctuation is 
clSp/p. In a statistically steady accretion disk, the local thermal time ~ {aQ)~^, so if the correla- 
tion time of mechanical fluctuations is ^ entropy fluctuations due to dissipation are expected 
to be at most oc a. Entropy fluctuations may also arise by advection of material in the presence 
of a background entropy gradient, 5S ~ — (5x • VS. The simulations cited above indicate that the 
turbulent Reynolds stress scales with the magnetic stress so 5v ~ a^/^c^, and we assume that the 
largest displacements have frequencies ~ Q, so the r.m.s. displacement (5x ~ a in any direction. 
Vertical and horizontal entropy gradients are likely to have scale lengths ~ h and ~ r, respectively, 
which would seem to make vertical displacements more important. However, the azimuthal force, 
and hence torque, exerted by a mass element changes only at second order as it is displaced verti- 
cally through 6z ~ a^^'^h <^ h, so vertical displacements contribute at 0(a). Radial displacements 
contribute at first order, (6Ti/T,)(is/dr ~ a^^'^h/r. The vertically uniform disk of N2005 has a radial 
entropy gradient, since the sound speed is constant but the density is not; however, since N2005 
obtains (a) 5 x 10~^ (in agreement with previous work), it happens that a^/^ 0.07 = h/r in 
these models, so that we cannot distinguish between ~ a and ~ a 

It would be useful to calculate torque fluctuations in a variety of disk models with varying 

entropy gradients, thicknesses, and (if possible) strengths of turbulence. Lacking such information, 
we provisionally assume 5T, ~ f{a)T,, where / is some function (probably linear). This leads to the 
paramctrization (7) if a is a universal constant. As will be seen, diffusion then dominates mean 
migration, i. e. inward drift at the rate predicted for a laminar disk, in the outer parts of plausible 
disk models. If ST, is also proportional to h/r, diffusion becomes even more dominant in the outer 
parts of a flaring disk. 

N2005 reports that the r.m.s. fluctuating torque per unit planetary mass is 1.5 ± 0.5 x 10~^ 
in code units. In the same units, 27rGSrp 3.3 x 10~^ (Nelson 2005, private communication); the 
latter is independent of radius because E oc r~^. Therefore, 



Cd = 0.046 



(8) 
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is taken as the reference value of the dimensionless factor in equation (7). 

Following Eq. (2), the diffusion coefficient D{J) is the product of the correlation time Tc 
(6) and the torque variance (7). Both factors are uncertain, but within the framework of this 
advection-diffusion approach, they combine into a single parameter. We hope that future numerical 
simulations will refine our estimates of this parameter and, more importantly, test the validity of 
the scalings (6) and (7). 

For comparison. Ward (1997a) 's expression for the mean Lindblad torque density, 
leads to the scaling 

For consistency with our calibration of torque fluctuations, wc also calibrate this relation against 
Nelson's simulations of laminar disks, in which the type I torque per unit mass is ~ 1.5 x lO~^Mp/M0 
in code units. The constant of proportionality in Eq. (10) is then 2± —4.8. Using Eq. (10) implies 
that the mean torque is always negative. This is fine for plausible power-law disks but might be 
misleading for alpha disks, where the surface density and temperature can be highly structured, 
allowing in principle for positive (or at least much reduced) T at some radii. Thus when modeling 
alpha disks we use Eq. (9) directly, as described in Menou &: Goodman (2004). 



3. Analytic results 

Power-law disk models such as the Minimum-Mass Solar Nebula (hereafter MMSN; Hayashi 
(1981)) allow the advection-diffusion equation to be solved analytically, at least in steady state. 
Some aspects of time-dependent evolution can also be determined analytically for this class of disk 
models. Throughout this section, all disk models are considered to be infinite in radius as well as 
steady. While unrealistic, this is technically convenient, and we are still able to identify results that 
are likely to be sensitive to the details of the outer boundary conditions: for example, the mean 
orbital lifetime. For reference, the MMSN model adopted here has mean surface density E cx r~^/^ 
and aspect ratio h/r cx r^l^. 



3.1. Steady state 

A steady-state distribution of planets, /(J), satisfies 

^ = s(A (11) 
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where the right-hand side is a source term representing the rate of formation of planets at radius rp 
with angular momentum J = Mpy^GM^r^. Analytic solutions to Eq. (11) exist when mean torque 
and the diffusion coefficient are both power-laws in J: 

f (J) oc D{J) oc J^, (12) 

and the source term is a Dirac delta function 

S{J) = 6{J - Js)JsA, (13) 

where A is the rate of planet formation at Js- We see immediately that the flux of planets interior 
(FJ) and exterior {Fj') to the source are both constant and related to each other by 

F+-Fj = JsA. (14) 

First we consider the limiting behavior as J ^ oo. The local timescale for secular inward migration 
is given by tmig = <^/|f (J)|, and the local diffusion timescale is given by tjifr = J"^ /D{J). The 
ratio, idiff/^mig oc J"-/3+i indicates that the mean torque becomes asymptotically negligible if 
a — 13 + 1 < Q. (This condition holds for the minimum-mass solar nebula, in which V oc —J"^ 
and D oc J.) For such disks the solution to Eq. (11) at large J exterior to the source becomes 
/(J) (C — JFj')/D{J), where C is a constant of integration. We insist that FJ > since 
there should not be an incoming flux of planets from infinity. Thus the numerator eventually 
becomes negative. The diffusion coefficient D(J) is unquestionably positive, so / becomes negative 
at sufficiently large J. But this is absurd, since f{J)dJ is a probability. Thus we cannot have 
both a positive flux and a positive / at large J in steady-state, and it follows that Fj = and 
FJ = — J5A. All the flux generated at the source is accreted through the inner edge of the disk. 



We now discuss Eq. (11) in full, emphasizing behavior at small and large J. Let 

dn __ f ( J) 

dJ ~ ~ D{J) ' 

so that exp[//(J)] is an integrating factor for Eq. (11). Then 



(15) 



(16) 



^D{Js)f{Js)-Ff J e^^'UJ . 

As before, the upper (lower) sign on and Ff applies to J > J5 (J < Js)- When V{J) and D{J) 
are power laws, can be evaluated explicitly, but we will not write it out. Normally the mean 
torque F < 0, so that d^/dJ > 0. Then we can choose the constant of integration in Eq. (15) so 
that < for J < J5 and > for J > Js- All disks of interest in this paper are dominated 
by the mean torque at small radii and by diffusion at large J. Hence exp[— /x(J)] becomes large at 
small J, and so for a well-behaved solution, the contents of the square brackets in Eq. (16) must 
tend to zero as J — >^ 0. This defines a relationship between /(J5) and FJ . For J > J5, e'^^''^ > 1, 
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so the integral increases at least as fast as J — J5. As before, if the disk is indefinitely extended 
then we must have Fj = to prevent f~^{J) from becoming negative at large J. 

The main conclusion of this discussion is that as long as the mean torque is nonpositive in the 
outer disk, the flux of planets diffusing to infinite radius vanishes. We have reached this conclusion 
by considering steady solutions but will assume that it applies generally. 



3.2. Time-dependent evolution 

The time-dependent problem is more difficult, but we can determine some general character- 
istics of the solution by analytic means, most importantly that the probability of surviving without 
accretion in an infinite steady disk, though declining monotonically with time, can have a long power 
law tail. For the special case of the MMSN, where D oa J and V oc —J~^, this probability declines 
only ast~^. The asymptotic form of the solution at late times is self-similar except near J = and 
is described explicitly in dimensionless units by equation (32). 

The evolution of a planet created with J = J5 at i = can be represented by solving Eq. (5) 
with a source term 5{t)5{J — Js). The Laplace transform 

00 

fis,J) = J e-''fit,J)dt, (17) 


reduces the time-dependent problem to 
52 



D{J)f{s, J)] - ^ [t{J)f{s, J)] - sf{s, J) = -S{J - Js). (18) 

Unfortunately, although Eq. (18) is just an ordinary differential equation in J, we cannot solve it 
completely even for power law disks described by Eq. (12). The difficulty is in part that physically 
interesting cases have a — P < —1, which is precisely the condition for Eq. (18) to have an irregular 
singular point at J = 0. However, a formal solution can be obtained in powers of s, 

00 / \ ^ 

/(^"^) = E^/-('^)- (19) 

n=0 

Replacing exp(— st) by its Taylor series in the Laplace transform (17) shows that fn{J) is the n*^ 
moment of f{t, J) with respect to t, 

OD 

fn{J)= J ef{t,J)dt. (20) 


More interesting than the moments of / are the moments of the corresponding flux: 

00 

F„(J) = -(D/„)' + f /„ = j t^Fjit, J) dt. (21) 
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The probability that the planet survives longer than t before accretion, P{t), is related to the flux 
at the origin by dP/dt = Fj{t,0). Therefore, 

oo 

r dP 

-Fn{0) = -J e—dt (22) 



is the n*'^ moment of the lifetime, and in particular, — i^i(O) is the mean lifetime. 

Thus we are motivated to find the fn{J)- Substituting Eq. (19) into Eq. (18) leads to 

[DiJ)MJ)]"-[fiJ)MJ)]' = -d{J-Js), (23) 
[DiJ)fn{J)]" -[fiJ)fn{J)]' = -nfn-i{J) n > 0, (24) 

where the primes indicate derivatives with respect to J. The solution to Eq. (23), which is formally 
identical to the steady-state case (11), is a Green's function for eq. (24). It can be constructed from 
homogeneous solutions yo{J) and yi{J) with fluxes and —1, respectively: 

[DiJ)yoiJ)]' -f(J)yo(J) = 0, (25) 
[D{J)y,{J)]' -f(J)yi(J) = 1. (26) 

The Green's function is 

G{J, Js) = yi{J<)yo{J>)/yo{Js), J< = min(J, Js), J> = max(J, Js). (27) 
The formal solution to the system (23)- (24) is then 
/o(J) = GiJ,Js), 

oo oo 

fniJ) = n\JdJi...J dJn GiJ, Ji)G(Ji, J2) . . . G(Jn, Js)- (28) 



The existence of the moments (20) therefore depends upon the convergence of the integrals in (28). 

An important special case is the MMSN, where D(J) = DqJ and T = — Tq It is convenient 
to choose = (Tq/Dq)^/'^ as the unit of angular momentum and = J^^/Dq as the unit of time 
so that D{J) J and f (J) —J~^. The homogeneous solutions (25) become 

2 1 

yo(J) = J-ig"^ where x = — ^ , (29) 

Jv 2 



yi 



(J) = 1 - 7r^/2xe-^'erfc(ar) 



r 72 



< 



J2-(3!!)J4 + (5!!)J6-... J^0+ 

00 ^ 00 ^ 

^ ~ \f\ ^ (2^)11'^ " + X) (2n-l)!!'^ " J ^ CO. 
n=0 ^ ' n=l ^ ' 



(30) 



The first expansion is only asymptotic, in keeping with the fact that J = is an irregular singular 
point of Eq. (26), whereas the second converges for all J 7^ 0. 
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Armed with these results, we can now determine the convergence of the integrals (28). Because 

yo{J) oc as J ^ cxd, it follows from (27) that the integrand for /i(J) is proportional to J^^ as 
Ji ^ oo. Thus the integral over Ji is logarithmically divergent. Consequently the mean lifetime 
defined by Eq. (22) is not finite. Since the divergence is only logarithmic, one may guess (correctly) 
that P{t) (X at late times in the infinite MMSN. 

The long tail in the distribution of lifetimes results from diffusion to large radii where torques 
are weak and planets linger. It is natural to seek a self-similar solution to Eq. (5) that reflects this 
behavior. Unfortunately J* defines a preferred scale that prevents self-similarity. Omitting T from 
Eq. (5) on the grounds that this term is unimportant at J S> J* removes the obstruction: 

%{t,J) = ^[Jf{t,J)], (31) 

in the present dimensionless units with D{J) oc J as in the MMSN. The tilde distinguishes the 
solution to the approximate equation (31) from the exact solution f{t,J) of Eq. (5). We require 
/ — > as J ^ oo. On the other hand, /(t, J) should be positive and finite at J = for the following 
reason. At late times such that J <C ^^D{J)t, we expect f{t, J) to match onto a multiple of the 
constant-flux steady-state solution defined by Eq. (30): f{t,J) ^ —Fj{t,0)yi{J). We may replace 
/ with / in this relation if J S> J*, and in that region yi{J) ~ 1. Eq. (31) can be solved with these 
boundary conditions by separation of variables.^ For arbitrary positive initial conditions, 

fit, J) oc t-^g-^A ast^oo. (32) 

As promised, this has self-similar form: J enters only in the combination J/t [or J /(Dot) in 
dimensionful units]. It can be verified that Eq. (32) is an exact solution of Eq. (31), and its integral 
over all J is oc t~^, in agreement with our expectation for P{t) at late times. 

For other values of the indices a k, (3 satisfying a < /? — 1, we have the following situation. 
Again, f is negligible at large J, so Eqs. (25) & (27) yield G{Ji,Js) oc l/D{Ji) at Ji » Js- 
Therefore the integral (28) for the first moment /i(J) diverges if /3 < 1, and consequently, the 
mean lifetime converges only if /3 > 1. In other words, the mean lifetime exists only if tdis increases 
more slowly than r^/^. 



3.2.1. A useful theorem about the eigenvalue spectrum 

The following technical result helps to interpret differences between the late-time behaviors of 
the planetary distribution in the MMSN and in alpha disks. 

The largc-t behavior of f{t, J) is dominated by the singularities of f{s, J) at largest Real(s). 
These correspond to eigenvalues {A} obtained by separating variables in Eq. (5) with assumed 



^This is eased by the change of variables x = \/2J and y{x) = xf{J), which produces the Bessel equation of order 
one on the righthand side. Hankel transforms can then be used. 
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time dependence exp(— At). If the eigenvalues are real and discrete, < Ai < A2 < . . ., then 
f{t, J) oc fx^{J) cxp(— Alt) at late times, where f\{J) is the cigcnfunction corresponding to A. (We 
assume that these modes arc complete.) On the other hand, if the eigenvalue spectrum is continuous 
and extends to arbitrarily small positive A, then the late-time decay is slower than exponential, 
perhaps a power law. Defining 

gx{J) = J fx{J)dJ, 
J 

and setting w = and p = e^D, with f^ as in Eq. (15), one can re-cast the separated version of 
Eq. (5) in self-adjoint Sturm-Liouville form 



_d_ 

dJ 



+ Xw{J)gx{J) = 0. 



(33) 



In all disks of interest to us, / ^ as J ^ and the integral of / over J is finite, whence 5^(0) 
and gx{oo) = 0. A theorem of Drabek &: Kufner (2005) then asserts that if the limit 



lim 

J—* 00 



w{x)dx 



1/2 



dx 



.1/2 



Uj 



p{x)_ 



(34) 



vanishes, then the spectrum is discrete and positive. [The theorem also requires p and w to be pos- 
itive and continuous, which is always true for us, and it encompasses "quasilinear" generalizations 
of Eq. (33).] It is easy to show that in disks where tmig/idiflf tends to zero as J — and to infinity 
as J — cxD, the limit (34) is equivalent to 



lim 

J— >oo 



J 



/; 



dJ 



.1/2 



(35) 



In the MMSN [D{J) oc J], this limit diverges, so a continuous spectrum extending to A = 0+ is 
allowed and indeed must exist, or else we would not have the self-similar behavior discussed above. 
In the alpha disks studied numerically in §4, on the other hand, it appears that tdiff = JV^(^) ^ 
as J ^ 00 [Fig. (1)], so the limit vanishes and a discrete spectrum is implied. This is clearly true 
of power-law disks with a — [3 + 1 < Q and /? > 2, so the survival probability P{t) should decay 
exponentially in these disks. We conjecture that the continuum exists and P{t) follows a power 
law for all /3 < 2, because further consideration of Eqs. (28) indicates that the n*^ moment of the 
lifetime exists only if /3 > 2n/(n -|- 1). 



4. Numerical method and tests 



The diffusion equation (5) with a source term reads: 



-gj=Sit,J)-^. (36) 
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We find numerical solutions with an implementation of the generic relaxation code described in 
Hameury et al. (1998). The finite-difference scheme is fully implicit in time. Given a distribution 
f{t, J) at time t, we only need to solve the one-dimensional boundary- value problem in J to find 
the solution at time t + dt. We thus relax successively from one solution to another over a variable 
timestep, dt. Controls over the size of the timestep ensure that, at each grid point, the change 
in each variable over dt is within an acceptable tolerance. The system is described by two first- 
order differential equations equivalent to Eq. (36), supplemented with a third differential equation 
that defines the numerical domain. For simplicity, we do not implement an adaptive mesh but 
use instead a fixed grid on the numerical domain, chosen to be uniformly spaced in logr. We 
take as variables the set {f,d{Df)/dJ, J). Numerical stability of the solutions is achieved with 
centered-differencing schemes for / and J, and an upwind-differencing scheme for d{Df)/dJ. 

Three boundary conditions are required. The numerical domain extends from the inner edge 
of the disk, at ~ 10~^ AU, out to 1000 AU. At the inner boundary we trivially set the value of J 
and are free to fix / = 0. At the outer boundary we demand the flux be zero. In steady-state, 
as discussed in §3.1, this is a boundary condition that must be satisfied if diffusion dominates at 
large radii. To avoid potential complications near the disk edge, we extend our numerical domain 
to sufficiently large radii that it does not influence our results. 

We have tested our numerical results for power-law models against the analytic steady-state 
solutions described in §3.1. The numerical solutions converge satisfactorily for a modest grid 
spacing (typically 1000 grid points over five decades in radius). Additional tests at ten times this 
resolution show clear numerical convergence of the steady-state results. Agreement with the late- 
time behavior expected from the time-dependent self-similar solutions described in §3.2 will also 
be demonstrated below (see Fig. 4). 

5. Numerical Solutions 

We study the combined effects of migration and diffusion on the orbital evolution of low-mass 
proto-planets for two separate classes of disk models: the Minimum Mass Solar Nebula (MMSN) 
and T Tauri alpha disks. The MMSN model adopted here has a mean surface density S oc r^^/^ and 
an aspect ratio h/r oc r^/^, which are normalized to 4200gcm~^ and 0.1 at r = lAU, respectively. 
The T Tauri alpha disk model adopted is identical to the irradiated disk described by Menou & 
Goodman (2004, see their Fig. 1), accreting steadily at a rate M = lO~^M0yr~^. We take as our 
reference model the alpha disk with viscosity parameter a = 0.02. For comparison, we also consider 
a more massive and cooler alpha disk, with a = 10~^. The inner edge of the disks is located at 
0.01 AU and 0.02 AU in the MMSN and the alpha disks, respectively. The mass of the primary is 
IMq in the MMSN and 0.5Mq in the alpha disks. These slight model differences have very little 
influence on our results. Note that, while the disk models used are assumed to be steady, we are 
studying the time-dependent orbital evolution of planets embedded in these disks with equation 
(5). As a consequence, our method can only apply to proto-planets of small enough mass that they 
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exert a negligible feedback on the disk structure (broadly speaking, Mp < lOM^; e.g. Menou & 
Goodman 2004). 

Turbulence is assumed to be present everywhere in our disk models. The magnitude of the 
corresponding diffusion coefficient is uncertain through both the amplitude and the correlation time 
of the turbulent torque fluctuations. We quantify these uncertainties with a single parameter, e. If 
we define 

then 

DiJ) = re X mt,J)f = e(27r)3(0.046)2^^^. (38) 

The calibration performed in §2.1 suggests that e ~ 0.5 is a reasonable value, which we have assumed 
to be independent of radius. However, we also investigate the effects of varying this parameter by 
several orders of magnitude. 



5.1. General Scalings 

It is useful to discuss a few general scalings specific to our disk models before we discuss the 

details of time-dependent numerical solutions. To quantify the relative importance of diffusion vs. 
migration in the models, we define a local diffusion time tjig = / D oc S~^J~^ oc J~4fc-5 
S oc r^) and a local migration time, tmig = J/|f | oc h?Ti^^ J^^ M^^ , according to Eq. (10), so that 
tmig oc J'^Mp^ in the MMSN. The ratio of these timescales has a particularly simple dependence 
upon disk and planetary properties: 

— 0^ t:F • (39) 
idis Mp 

One expects the evolution of planets located in regions of a disk with tmig <^ tdiff to be 
dominated by migration and vice-versa. While the diffusion time is independent of planet mass, 
the migration time is inversely proportional to Mp. Therefore, lower mass proto-planets will be 
more sensitive to the effects of turbulent torque fiuctuations. In addition, the two timescales differ in 
their dependence on the mean surface density, S, which tends to make planets in denser disks more 
susceptible to the effects of diffusion. As we shall now see, a general consequence of these scalings 
is that one expects the effects of diffusion to dominate only in the outer regions of proto-planetary 
disks (at least for the classes of disk models considered here). 

Figure 1 shows profiles of various important quantities in the MMSN and our two alpha disk 
models. Panel I shows profiles of mean surface density, S. Both the slope and normalization of 

E profiles matter for the relative importance of migration vs. diffusion. Panels III and IV show 
profiles of local migration and diffusion times in the various disk models. While tjifr increases 
linearly with J (oc r^^^) in the MMSN, it decreases with increasing radius in the alpha disk models. 
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This is a consequence of the steeper S profile in the MMSN. The profiles of tmig, derived from 
Eq. (10) for the MMSN and from the detailed calculations'^ of Mcnou & Goodman (2004) for the 
alpha disk models, show local migration times which decrease toward small radii in both classes of 
models. Finally, panel II shows the ratio of timescales, tmig/^diff > which is important to characterize 
the behavior of planets embedded in these disks. As mentioned earlier, in all three models, this 
ratio grows above unity (indicated by the horizontal dash-dotted line in panel II) only at radii 
> 10 AU if Mp < 1 M0 , indicating that diffusion should dominate over migration only in the outer 
regions of these disks (and this despite the increasing value of tdiff with radius in the MMSN) . 

Although very useful, the profiles shown in Figure 1 do not capture all of the complexity of the 
migration-diffusion problem. First, all the quantities plotted in panels II-III-IV were derived for 
a planet mass Mp = IM^ and a torque fluctuation normalization e = 0.5. The profiles are easily 
rescaled for different values of these parameters (as indicated by labels in Fig. 1), but it is clear 
that the exact location of the transition from migration- to diffusion-dominated evolution depends 
on the proto-planetary disk model adopted, the value of the planet mass, and the normalization 
e of the diffusion coefficient. Nevertheless, the dominance of diffusion over migration in the outer 
regions of these disks is robust because S/i^ tends to increase with radius in disk models of interest 
[Eq. (39)]. 

Secondly, as shown in panel III, in a disk model with realistic opacities, sudden and non- 
monotonic variations in the value of tmig at radii < 3 AU could complicate the orbital evolution of 
planets substantially, by stalling their migration as they reach the disk inner regions. This is best 
explored with global numerical solutions of the time-dependent Fokker-Planck equation, as seen in 
§5.2. 

A third complication arises from the radial dependence of the diffusion driven by turbulent 
torque fluctuations. As mentioned in §2, the radial variation of the diffusion coefficient leads to 
an additional contribution to the advected ffux of planets. Since —dD/dJ oc —{7 + 4k) J^'^'^'^ 
for E oc r^, and since k = —3/2 in the MMSN and k > —3/2 in the alpha disk models, this 
additional term tends to advect planets inward in all the disk models considered. We can define a 
local advection timescale associated with this new contribution, J/{dD/dJ), and compare it to the 
values of tmig and idiff defined above. We find that this advection timescale is everywhere exactly 
equal to t^is in the MMSN, while it is shorter than tdis everywhere in our alpha disk models (as 
shown explicitly for the a = 0.02 case by the short-dashed line in panel IV of Fig. 1). Consequently, 
even in disk regions where diffusion nominally dominates over migration in our models (as measured 
by the ratio ^diff/^mig), this additional contribution due to radially inhomogeneous diffusion should 
still cause a rather significant effective advection of planets towards the star. 



^Specifically, we use the 3D model in which the gas is assumed to be vertically extended over a full disk scale 
height for the torque calculation. 
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5.2. Continuous and Burst Formation of Planets 

We have studied two different types of time-dependent models for the diffusive migration of 
planets in proto-planetary disks, corresponding approximately to a continuous and a burst scenario 
for planet formation. In both cases, for the sake of clarity, the source of planets was localized, 
taking the form of a sharply peaked Gaussian approximating the delta source function of Eq. (13). 
We have verified that our results do not depend sensitively on the detailed shape of the Gaussian 
source function adopted. As we shall see, the predictions for these rather different scenarios for 
planet formation are broadly consistent with each other. 

Figure 2 shows steady-state distributions obtained by numerically solving Eq. (36) with a source 
term centered at rg = 10 AU, steadily producing earth-mass planets at a rate A = (10^ yr)"^- The 
topmost curve in each panel is the steady-state distribution. The enclosed curves, on the other 
hand, are snapshots of the decaying distributions at selected time intervals after the source is shut 
off. We have verified that the steady-state distributions for the MMSN in Panels I and III are fully 
consistent with the analytic solution, Eqs. (15)-(16), for power-law disks. A fiducial normalization 
of e = 0.5 for torque fluctuations is assumed in all cases except in panel III, where e = 0.02 has 
been scaled down to extend the migration-dominated regime out to ~ 10 AU in the MMSN. Thus, 
at 10 AU, the disks in Panels I and IV are diffusion-dominated and the disks in Panels II and III 
are migration-dominated. 

The distinction between diffusion-dominated and migration-dominated evolutions is most clearly 
illustrated by comparing panels I and III. While the initial steady-state peaked distribution is eroded 
and advected for the most part in panel III, with little diffusion at large radii, there is significantly 
more diffusion at large radii in panel I. That distinction is not as clear in the alpha disk models, 
however. Note that, in both panels II and IV, the initial steady-state distribution is strongly peaked 
not at the source radius, but at radii where migration stalls (as seen from panel III in Fig. 1). In 
addition, even though diffusion is more important in the model with a = 10~^ in panel IV, the 
effect appears to be minor in smoothing out only slightly the same sharp distribution features as in 
panel II. This is partly explained by the strong contribution of radially inhomogcneous diffusion to 
inward advcction of planets, which is always important in our alpha disk models (as explained in 
the previous subsection). Even though diffusion is equally important at rs = 10 AU in the models 
of panel I and IV (see tmig/idiflf in Fig- 1)) it is more inhomogcneous in the alpha disk model and 
thus promotes a significantly more efficient inward advection of planets. 

Our models with a burst of planet formation confirm these general trends. To study this 
formation scenario, we solve the diffusion equation (5) without source term and impose a specific 
initial condition for the distribution of planets at t = instead. The initial distribution, f{t = 
0, J) = fi{J), is modeled as a narrow Gaussian centered on Js = Mp\/GM^rs- Figure 3 shows the 
evolution of such distributions of Earth-mass planets initially located at rg = 10 AU in the same 
four disk models as in Fig. 2. 

As before, a comparison between the two MMSN models shows that the effects of diffusion are 
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much more important in the model with e = 0.5 (panel I) than in the one with e = 0.02 (panel III). 
In particular, a sub-population of planets diffusing at large radii is clearly visible in panel I. On 
the contrary, while diffusion still acts to widen the initial distribution in panel III, the dominant 
effect is a global advection of the distribution of planets. Finally, in agreement with our previous 
results, Figure 3 shows that the most significant feature of alpha disk models is the effect of rapid 
variations in the migration rates of planets, which causes them to stall at specific locations in 
the disk and results in time-varying peaks in the distributions of panels II and IV. Even though 
diffusion is comparatively as important in panels I and IV, at 10 AU, its effects appear to be much 
less important in the alpha disk model of panel IV. We attribute this difference, once again, to the 
more inhomogeneous diffusion in alpha disk models, which contributes substantially to the inward 
advection of planets into regions of the disk where diffusion becomes gradually less important and 
the effects of stalling radii are felt more strongly. 



5.3. Radial Excursions and Lifetimes 

Two important aspects of the orbital evolution of planets which are not well captured by 
Figure 3 are the properties of the surviving minority of planets at late times and large radii, and 
the representative planetary lifetimes in burst-formation models. 

The fraction of planets that diffuse to orbital radii greater than their formation radius, rs, can 
be characterized by the time-dependent quantity 

OO / oo 

Xa{t) = J f{t,J)dJ / J f{0,J)dJ, (40) 

where Ja is the angular momentum at a fiducial radius a[Au] > rs- Hereafter, the integral over 
the initial distribution appearing in the denominator of Eq. (40) will be normalized to unity. The 
maximum values of X2o{t) and XiQ^it)-, reached typically just after the majority of planets were 
accreted by the star, are listed in Table 1 for the four models of Figure 3. Not surprisingly, models 
with a strong role for diffusion tend to have a fairly large maximum fraction of planets diffusing 
beyond 20 AU 10%) and 100 AU 1%), while models with little diffusion have much reduced 
maximum fractions of planets at large radii (< 0.01% beyond 20 AU, < 10~® beyond 100 Au). 

Table 1: Maximum fraction of planets beyond a = 20 AU and 100 AU in the four models of Figure 3 



Disk Model 




Torque Fluctuations 


A20 




, ,max 
AlOO 




MMSN 




e = 0.5 


1.4 X 10" 


-1 


3.4 X 10" 


-2 


Alpha-Disk (a 


= 0.02) 


e = 0.5 


1.9 X 10" 


-4 


1.4 X 10" 


-6 


MMSN 




e = 0.02 


1.0 X 10" 


-4 


7.6 X 10" 


-8 


Alpha-Disk (a 


= 10-3) 


e = 0.5 


7.4 X 10" 


-2 


9.7 X 10" 


-3 



We compute three additional time-dependent quantities to help us better characterize the 
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properties of surviving planets in our models. The first is the survival probability, computed as 

oo 

P{t) = J fit, J)dJ. 



The second is the median orbital radius of surviving planets, rmedC*)) defined by 

M"p-y/GM,rmed(t) 

f{t,J)dJ=^-P{t). 


Finally, since the mean lifetime is often ill-defined (§3.2), we calculate the median lifetime tm defined 
by P{tm) = 1/2. 

Figure 4 shows how the survival probability and the median radius of surviving planets evolve 
with time in burst formation models for Earth-mass planets in the MMSN. Three separate models 
with burst formation at rs = 1, 10 and 100 AU are shown for comparison (solid lines). Note that 
the middle solid line corresponds to the exact same model as the one shown in panel I of Figure 3. 
For reference, vertical dotted lines bracket typical lifetimes for gaseous proto-planetary disks (10®- 
10^ yrs; e.g. Haisch et al. 2001). The median lifetime of the population of planets, which is such that 
P{t) = 0.5, is indicated by a horizontal dash-dotted line in the upper panel. The lifetimes differ 
substantially in the three models shown, from much shorter to somewhat longer than the typical 
disk lifetimes. The self-similar P{t) oc behavior predicted in §3.2 is indeed achieved at late 
times in the three models. As shown in the lower panel of Figure 4, the median radius of surviving 
planets does not evolve much at early times and remains close to the formation radius, as the 
distribution of planets mostly diffuses out to both small and large radii. Once a significant fraction 
of planets have been accreted by the central star, however, the surviving planets are found at a 
median radius increasing almost quadratically with time independently of where they were formed. 
This is just what one would expect from the self-similar solution (32): r^^^ oc Jmed Dgt ln2. 

Figure 5 shows the evolution of the same two quantities in burst formation models for Earth- 
mass planets in the alpha disk (with a = 0.02) this time. Six separate models with burst formation 
at rs = 1, 5, 10, 20, 100 and 500 AU are shown for comparison (solid lines). Note that the fourth solid 
line from the top corresponds to the exact same model as the one shown in panel II of Figure 3. 
The results are qualitatively consistent with the MMSN ones, with a few important differences. 
There is no self-similar behavior at late times in this case and the more dominant role of advection 
in the alpha disk model tends to make the median radius of surviving planets move inward initially, 
until the majority of them are lost to the central star. Among the surviving minority, however, 
rjncd asymptotes to ~ 50 AU at late times, independently of the formation radius, rs- This behavior 
is explicable if, as predicted by the theorem cited in §3.2.1, the eigenmodes of Eq. (5) are discrete: 
then f{t, J) converges upon the most slowly decaying mode, and rmed(^) tends to the median radius 
of this eigenmode. Figures 1 and 5 suggest that this radius is close to that at which -|- 
[more accurately, + f^^, where tgdv = J/{dD/dJ)] is minimized, as might be expected since 



/ 



-20- 



planets should linger longest where the combined effects of diffusion and migration are weakest. In 
the MMSN, both timescales increase without limit as J — oo, and the spectrum of eigenmodes is 
apparently continuous. 

It is also worth emphasizing that the largest peak value of fmig at sub-AU distances in alpha 
disk models (see panel III in Fig. 1) effectively limits the migration time of Earth-mass planets into 
their star to values > 10^ yrs (for a = 0.02), which exceeds typical disk lifetimes. 

Based on this general understanding of diffusive migration in our burst-formation models, let 
us now focus more systematically on several important model dependencies. Figure 6 illustrates 
how the median lifetimes of Earth-mass planets formed at different radii in the MMSN and the 
alpha disk are affected by changes in the normalization of torque fluctuations, e. For the nominal 
e = 0.5 value adopted (indicated by a vertical dash-dotted line in Fig. 6), diffusion mostly affects 
the orbital evolution of planets formed at the largest radii, by reducing somewhat their median 
lifetime. As e is increased above this nominal value, the trend becomes clearer. Median lifetimes 
for all formation radii are considerably reduced at larger values of e, more so in the MMSN than 
in the alpha disk. The long-dashed line in Figure 6 results when we neglect the mean torque F in 
the governing equation (5). The convergence of the solid curves for the alpha-disk model to a value 
independent of rs is related to the value of the diffusion time, tdiffj at the disk inner boundary. For 
sufficiently large e, tdiff ^ imig everywhere on the grid. But tdiff increases toward small radius in 
the alpha disks (Fig. 1). So the median lifetime in the large-e regime is comparable to tjifj at the 
inner boundary, Tmin; we have verified that tmed decreases as rmin is increased. Nevertheless, in all 
cases, the effect of diffusion is to reduce the median lifetime. 

Figure 7 illustrates how the median lifetimes of planets of various masses, all formed at rs = 
10 AU in the MMSN and in the alpha disk, are affected by changes in the normalization of torque 
fluctuations, e. For the nominal e = 0.5 value, diffusion mostly affects the orbital evolution of the 
lower mass planets, and reduces somewhat their median lifetimes. This results from the larger 
values of tmig at lower planet mass, whereas tjifj is independent of planet mass. As e is increased, 
diffusion progressively dominates and the distinction among planets of various masses gradually 
disappears. For large enough values of e, the median lifetime of all planets with masses My < lOM^ 
formed at rs = 10 AU is reduced below typical disk lifetimes. As in the previous figure, however, 
the lifetimes for e S> 1 in the alpha disk are dominated by the diffusion time at the inner boundary. 

6. Discussion and Conclusion 

As anticipated by Nelson 8z Papaloizou (2004), Laughlin et al. (2004), and Nelson (2005), 
our Fokker-Planck treatment confirms that diffusion driven by turbulent torque fiuctuations in 
proto-planetary disks can greatly influence the orbital evolution of embedded low-mass proto- 
planets. Contrary to some expectations, diffusion does not promote planetary survival but instead 
systematically reduces the lifetime of most planets in the disk. However, it does help a small 
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fraction of planets diffuse to large radii where they can survive for extended periods. (We consider 
that a planet "survives" if its orbit remains larger than the inner edge of the disk.) 

Our results point to a potentially important role for non-deterministic effects in planet forma- 
tion scenarios. In models with a significant level of diffusion, most proto-planets in the Earth-mass 
range could end up being accreted by their stars, while only a small fraction of them would survive 
by diffusing to large radii. In scenarios where all proto-planetary disks are efficient at forming 
planets, this could be interpreted as leading to only a fraction of all potential host stars being 
successful in keeping a system of surviving planets. Even systems starting with very similar initial 
conditions may end up with very different planetary orbital configurations once their gaseous disks 
disappear. 

Neglecting the effects of direct gravitational interactions between proto-planets, which are 
likely to happen from differential rates of migration and diffusion, is one of many model limitations 
in our work. Wc have already mentioned that our study is restricted to proto-planets of low enough 
masses that they do not affect in any way the structure of their host disk. This excludes the late 
stages of giant planet formation, which are obviously of considerable interest. In particular, it 
has been suggested that migration and diffusion, each independently, could accelerate the rate of 
growth of cores until they reach the critical mass at which runaway envelope accretion occurs (e.g. 
Rice & Armitage 2003; Alibert et al. 2005a,b). Studying these possibilities within the context of 
global models such as ours would be interesting. 

Even within the strict regime of applicability of our models, significant uncertainties exist. Let 
us mention a few. Wc have focused our work on strictly circular orbits. We have assumed that the 
underlying proto-planetary disk structure does not evolve with time, and that the amplitude 
of surface-density fluctuations is constant with radius. We have assumed that torque fluctuations 
can simply be added to a laminar mean torque which is not modified by the turbulence. The 
laminar torque itself is calculated from a simple extension of the two-dimensional Lindblad torque 
theory, which does not account for the role of corotation torques. 

Despite all these uncertainties, we believe that our results should be broadly valid, provided 
that turbulence and torque fluctuations in proto-planetary disks satisfy the various stochasticity 
criteria enunciated in §2. Understanding better the nature of turbulence in proto-planetary disks 
and in particular indications for long-term correlations in MHD simulations may turn out to be 
crucial in this context. We would like to issue a plea to the simulation community to study these 
questions, and in particular, whether the dimensionless amplitude and correlation time of density 
and torque fluctuations depend only upon the strength of the turbulence {via a) or upon the 
thickness h/r as well (§2.1). The answers will matter both quantitatively (in determining the 
relative importance of migration vs. diffusion in global models such as ours) and qualitatively (in 
justifying a Fokker-Planck approach or in falsifying it). 

We have seen that details of the protoplanetary disk structure (MMSN vs. alpha disks) do 
influence our results at the quantitative level. In that respect, the possibility that a magnetically 
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layered structure with quiescent dead zones exists in weakly-ionized regions of proto-planetary 
disks (Gammie 1996; Glassgold ct al. 1997; Balbus & Hawley 1998; Fromang et al. 2002; Ilgner 
&: Nelson 2006) may be important. Indeed, the results of Fleming &: Stone (2003) suggest that 
the nature of turbulence in dead zones differs qualitatively from that in MHD turbulent regions, 
perhaps leading to comparatively reduced torque fluctuations (see also Matsumura & Pudritz 2003, 
for other suggested effects of dead zones in planet formation scenarios). 

The general tendency for a fraction of planets to diffuse to large radii in proto-planetary disks 
may have consequences of some observational relevance. We note, for instance, that Veras & 
Armitage (2004) have looked into scenarios for outward planet migration in an attempt to explain 
observed structures in dust disks as resulting from perturbations by planets located at radii beyond 
those at which they are thought to be able to form. Rather than invoking outward migration, 
this may be expected in scenarios with a prominent role for diffusion, as illustrated by our various 
models. 
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Fig. 1. — In all four panels, quantities were plotted assuming a planet mass Mp = IM® and torque 
fluctuations normalized to e = 0.5. The long-dashed lines correspond to the MMSN, while the solid 
and dotted lines correspond to the a = 0.02 and a = 10~^ T Tauri alpha-disk models, respectively. 
Panel I shows mean surface density profiles. Panel II shows ratio of local migration time to local 
diffusion time in the three models. The horizontal dash-dotted line indicates the formal transition 
between migration-dominated and diffusion-dominated regime for Earth-mass planets. Panels III 
and IV show local migration times and local diffusion times, respectively. The short-dashed line in 
Panel IV corresponds to the local advection timescale due to inhomogeneous diffusion, J/ [dD / dJ) , 
for the a = 0.02 alpha-disk model. 
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Fig. 2. — Snapshots of the decay of steady-state distributions after souree shut-off. In each case, 
the source formed Earth-mass planets at rs = 10 AU at a rate A = (lO^yr)"^ until shut-off. The 
disk model in Panels I and III is the MMSN with torque fluctuations normalized to e = 0.5 and 
0.02, respectively. The disk models in Panels II and IV are alpha disks with viscosity parameters 
a = 0.02 and 10~^, respectively (e = 0.5). Snapshots are shown at the following times after source 
shut-off, in log-years: (I) 5,5.5,6,6.5,7,7.5; (II) 6.8,7,7.2,7.4; (III) 5.8,6,6.2,6.4; (IV) 5.8,6,6.2,6.4. 
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Fig. 3. — Snapshots of the evolution of initiahy spiked distributions corresponding to a burst of 
Earth-mass planet formation at rs = 10 AU. The disk model in Panels I and III is the MMSN with 
torque fluctuations normalized to e = 0.5 and 0.02, respectively. The disk models in Panels II and 
IV are alpha disks with viscosity parameters a = 0.02 and 10~^, respectively (e = 0.5). Snapshots 
are shown at the following times after burst formation, in log-years: (I) 4,5,6; (II) 6.8,7,7.2,7.4; 
(III) 5.8,6,6.2,6.4; (IV) 4,5,5.4,5.8,6.2,6.4,6.5,6.6. 
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Fig. 4. — Evolution of the survival probability, P (upper panel), and of the median radius of 
surviving planets (lower panel) for Earth-mass planets formed in a localized burst at t = in 
the MMSN. Torque fluctuations were normalized to e = 0.5. In each panel, the three solid lines 
correspond to planets formed at rs = 100, 10 and lAU (from top to bottom). Vertical dotted 
lines bracket typical protoplanetary disk lifetimes. The horizontal dash-dotted line in the upper 
panel corresponds to P{t) = 0.5, which defines the median lifetime of planets. A declining fraction 
of planets survives for long times at increasingly large radii, in agreement with the self-similar 
P{t) (xt~^ scaling predicted analytically. 



-29- 



o 




-10 



-15 



10'^ 108 
t (yr) 



10 



10 



102 ^ 



7] 
• I — I 



CD 

6 



10-2 





1 1 1 M M 


1 1 1 1 M M 


1 1 1 1 M M| 1 


1 1 1 M M| 1 


- 
Illttt 










































































1 





1 





1 1 


1 1 






10' 



10^ 



10"^ 108 
t (yr) 



10^ 



10 



10 



Fig. 5. — Evolution of the survival probability, P (upper panel), and of the median radius of 
surviving planets (lower panel) for Earth-mass planets formed in a localized burst at t = in the 
T-Tauri alpha-disk model (with a = 0.02). Torque fluctuations were normalized to e = 0.5. In each 
panel, the six solid lines correspond to planets formed at rs = 500, 100, 20, 10, 5 and 1 AU (from 
top to bottom). Vertical dotted lines bracket typical protoplanetary disk lifetimes. As before, a 
declining fraction of planets survives for long times but in these models the median radii reach 
an asymptotic value at which the sum of the secular and diffusive contributions to advcction is 
minimized. The advection of most planets at intermediate times is clearly visible in the median 
radius evolution. 
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Fig. 6. — Median lifetimes as a function of torque fluctuation normalization, e, for Earth-mass 
planets formed in a burst at radii rs = 100, 20, 10, 5 and 1 AU (from top to bottom) in the MMSN 
(dashed lines) and in the alpha-disk model (solid lines; a = 0.02). The vertical dash-dotted line 
corresponds to our fiducial calibration of torque fiuctuations (e = 0.5). Horizontal dotted lines 
bracket typical protoplanetary disk lifetimes. The long-dashed line indicates the median lifetimes 
expected in the alpha-disk model in the limit of negligible mean torque. Diffusion is unimportant 
at low e values, but it can significantly reduce the median lifetimes of planets formed at any radius 
for large enough e values. 
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Fig. 7. — Median lifetimes as a function of torque fluctuation normalization, e, for planets of various 
masses all formed in a burst at rs = 10 AU in the MMSN (dashed lines) and in the alpha-disk model 
(solid lines; a = 0.02). The various curves correspond to planets of mass Mp = 0.01, 0.1, 1 and 

lOAf0 (from top to bottom). The vertical dash-dotted line corresponds to our fiducial calibration of 
torque fluctuations (e = 0.5). Horizontal dotted lines bracket typical protoplanetary disk lifetimes. 
Diffusion is unimportant at low e values, but it can significantly reduce the median lifetimes of 
planets of any mass (< lOM®) for large enough e values. 



